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MAXIMUM MARGINAL LIKELIHOOD ESTIMATION OF A MONOTONIC POLYNOMIAL 
GENERALIZED PARTIAL CREDIT MODEL WITH APPLICATIONS TO MULTIPLE GROUP ANALYSIS 


Abstract 


We present a semi-parametric approach to estimating item response functions (IRF) use- 
ful when the true IRF does not strictly follow commonly used functions. Our approach 
replaces the linear predictor of the generalized partial credit model with a monotonic 
polynomial. The model includes the regular generalized partial credit model at the low- 
est order polynomial. Our approach extends Liang’s (2007) method for dichotomous 
item responses to the case of polytomous data. Furthermore, item parameter estimation 
is implemented with maximum marginal likelihood using the Bock-Aitkin EM algo- 
rithm, thereby facilitating multiple-group analyses useful in operational settings. Our 
approach is demonstrated on both educational and psychological data. We present sim- 
ulation results comparing our approach to more standard IRF estimation approaches 
and other non-parametric and semi-parametric alternatives. 


Keywords: Item response theory, semi-parametric models, monotonic polynomial, 
item response function 
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1 Introduction 

Under many unidimensional item response theory (IRT) applications, it is common 
to assume that the underlying latent variable, 0, is normally distributed, and that item 
response functions (IRFs) follow one of many standard functions (e.g., Thissen & Stein- 
berg, 1986). For example, it may be hypothesized that the relationship between the latent 
trait and probability of item response to dichotomous items has a logistic shape (i.e., 2- 
parameter logistic or 2PL), which provides a close approximation to the normal ogive 
model (Birnbaum, 1968). IRFs for ordered polytomous items may use the binary logistic 
function as “building blocks” (Thissen & Steinberg, 1986) to construct the graded re- 
sponse model (Samejima, 1969) or generalized partial credit model (GPC; Muraki, 1992). 

In some cases, the assumption of a normally distributed latent trait or the assumption 
that the IRF follows one of these standard functions is violated. Proceeding with a 
standard approach may result in poor recovery of the true IRF and/or underlying latent 
trait estimates (e.g., Ramsay & Abrahamowicz, 1989; Liang, 2007). Methods have been 
proposed to estimate IRT models in the presence of each of these assumption violations, 
with it possible to model either a non-normally distributed latent trait or non-standard 
IRF, but not both simultaneously. For instance, if all non-standard IRFs tend to follow 
the same shape, researchers may suspect that a non-normally distributed latent trait is 
the culprit and choose to model the distribution using an empirical histogram approach 
(Mislevy, 1984; Woods, 2007a), Ramsay curves (Woods & Thissen, 2006; Woods, 2006, 
2007b, 2008) or Davidian curves (Woods & Lin, 2008). Alternatively, non-standard IRFs 
can be modeled using techniques such as non-parametric methods (Ramsay, 1991; Rossi, 
Wang, & Ramsay, 2002; Samejima, 1977, 1979, 1984; Sijtsma, Debets, & Molenaar, 1990), 
semi-parametric methods (Liang, 2007; Ramsay & Winsberg, 1991), or by using Bayesian 
non-parametric estimation (Duncan & MacEachern, 2008, 2013; Miyazaki & Hoshino, 
2009; Qin, 1998). 


The purpose of this research is to present a semi-parametric technique for model- 


October 30, 2014 és) 


ing non-standard IRFs that may have potential advantages over some existing methods. 
This work builds upon that of Liang (2007), who used a logistic function of a monotonic 
polynomial (MP) to estimate non-standard IRFs for dichotomous items. She developed 
an estimation technique that used provisional estimates of @ for each individual’s latent 
trait (similar to Ramsay, 1991, 2000) in the complete-data log-likelihood. In this work, we 
present a multivariate logistic function of the MP that can model ordered polytomous 
items, which we refer to as the GPC-MP model since it reduces to the GPC model at 
the lowest order polynomial. When the number of categories is two, the model reduces 
further to the logistic function of a monotonic polynomial. When both of these condi- 
tions occur, the model reduces to the familiar 2PL. Furthermore, we utilize Bock-Aitkin 
(Bock & Aitkin, 1981) maximum marginal likelihood estimation using the expectation 
maximization algorithm (EM MML). 

Use of the MP ensures that IRFs (or for polytomous data the basic building block 
functions) are monotonically increasing, which is not necessarily the case under other 
non-parametric models (e.g., Ramsay, 1991; Rossi et al., 2002).!_ Our approach does 
not require provisional @ estimates or information from an external test on the subject’s 
ability (e.g., Liang, 2007; Ramsay, 1991; Samejima, 1977, 1979, 1984). Furthermore, use 
of EM MML facilitates mixing of different item types such as MP items and other item 
types within the same test, as well as multiple group analyses in which group means are 
compared or differential item functioning is investigated - features that are not always 
readily available under other approaches. Finally, to our knowledge this is the first 
approach to allow semi-parametric estimation of IRFs for polytomous items. Although 
several non-parametric techniques are available for polytomous items (Abrahamowicz 
& Ramsay, 1992; Mazza, Punzo, & McGuire, 2013; Ramsay, 2000; Santor, Ramsay, & 
Zuroff, 1994; Santor, Zuroff, Ramsay, Cervantes, & Palacios, 1995; van der Ark, 2007), 


'We assume that monotonicity is desirable in many testing situations where a correct item will always 
indicate higher (or equal) ability for all regions of the latent trait. However, we do note that releasing 
constraints on monotonicity may be useful for probing for severe departures from monotonicity or when 
non-monotonicity is actually predicted. 
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the performance of such approaches are typically not compared to each other and, to 
our knowledge, do not include formal tests for differential item functioning. 

The remainder of this manuscript is therefore organized as follows. Section 2 de- 
scribes our parameterization of the MP and the GPC-MP for estimation of IRFs. Section 
3 describes the details of our estimation procedure. In Section 4 we illustrate three poten- 
tial applications of the GPC-MP model on real-world datasets. In Section 5 we present 
two simulation studies that examine test conditions under which the GPC-MP model 
performs better at recovering the true IRF versus standard approaches and a compar- 
ison of the GPC-MP model using EM MML against the alternatives of using Liang’s 
(2007) estimation procedure or a non-parametric kernel smoothing approach (Ramsay, 
2000; Mazza et al., 2013). Finally, Section 6 provides concluding remarks. 

2 The Proposed Item Response Model 
2.1 Monotonic Polynomial 
The basic building block for the model we present is a monotonic polynomial and its 


first derivative of the form (e.g., see Liang, 2007; Heinzmann, 2005, 2008): 


m(6|¢,b) = b10 b>67 Jp... 4 bb (1) 


m'(6|a) = ag + 410 + 4207 +--+ + ap,67* (2) 


Thus, ¢ is an intercept parameter, and b = [by,...,b2444] and a = [ao,..., ax] are coeffi- 
cients with b; = a;_1/t fort = 1,...,2k +1. The degree of the polynomial is controlled 
by a user specified value for k € [0,0o). For example, k = 2 would represent a 5th 
order polynomial. In our application, @ represents the latent trait. In one of its initial 
applications, the monotonic polynomial was used to estimate an unknown cumulative 
distribution function (CDF), with the monotonicity of the polynomial ensuring that the 
CDF be monotonically increasing (Heinzmann, 2005, 2008). The monotonicity require- 


ment was established by ensuring that the polynomial was of an odd-order (by using 
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2k +1), and had a positive derivative m'(6|a) for all 6. Elphinstone (1985) showed that 
the latter requirement could be accomplished by re-parameterizing the coefficients a in 


terms of A, « = [a 1,...,a%], and B = [B1,..., Bx]: 


AT Tai (1 — 20n0 + (07, + Bn)0*) ifk > 0 
m'(6|A, 0, B) = Tas ( nO + (a7, + Bn)O) 3) 
A ifk =0 
and implementing the constraints that all 6B > 0 and A > 0. To later allow unconstrained 


estimation, we re-parameterize /1,...,6, = exp(T),...,exp(™) and A = exp(w), re- 


sulting in: 
; exp(w) [Tf (1 — 2a, + (a7 +exp(t,))62) if k > 0 
m'(0|w,«,T) = (4) 
exp(w) ifk =0 
Liang (2007) noted that the coefficients for m'(6|A,«,B), a, = [ag +++ ao, |’, fora 


polynomial of degree 2k + 1 could be represented in recursive form as ax = Tyay_1 = 
Ty Ty_1 °° T2T1A with any given matrix T, containing parameters only corresponding 
to a, and Bx, and a special case being that ag = A. Under our re-parameterization, each 
matrix T; contains only parameters a, and % and the expression for the coefficients a, 


becomes: 


ay = Tpag_1 = TyTy_1 -++ ToT, exp(w) (5) 
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The matrices T;, have dimensions (2k + 1) x (2k — 1) and are of the form: 


1 0 0 0 0 0 
— 2a | 0 0 0 0 
az + exp(T%) —2a, 1 0 0 0 
oe 0 a2 +exp(%|%) —2az 0 0 0 
0 0 O +++ ae +exp(T) — 2x, 1 
0 0 QO «+: 0 a2 + exp(T) —2ny 
0 0 O «+ 0 0 a? + exp(T) 


The advantage of the recursive form for a;, as described by Liang (2007), is that it is 
relatively easy to code a computer program to construct the polynomial m(6|€,w,«,T) 
and its derivative m’(0@|w,a,t) on the fly up to an arbitrary order. We start at k = 0 
to compute ag. We then move on to k = 1 by a; = Tjao, and so on, working our way 
up until the desired k. The coefficients b,; necessary for constructing the polynomial 
can then be computed from ax. If any derivatives of a, with respect to the parameters 


tT and w are required, we only need to differentiate and save the corresponding matrix 


dag _ 
’ Oa, 


T that contains such parameters. For instance T3725! exp(w), since «1 is only 
contained in Tj. 
2.2 A Monotonic Polynomial Generalized Partial Credit Model 

Suppose we have i = 1,2,...,N independently sampled individuals respond to j = 
1,2,...,n items, and c = 0,1,...C; — 1 coding the response categories for item j. When 
items are dichotomous (i.e., C; = 2), the two-parameter logistic (2PL) model can be used 


to represent the probability that individual i will endorse category 1 on item j: 


1 
4a exp(—(d; + 1;9i)) 
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where 6; is the person’s latent trait, 7; is a slope parameter, and 6; is the intercept. 
With polytomous items (C; > 2), the generalized partial credit model can be used to 
represent the probability that respondent i’s response to item j is category c (Muraki, 


1992), and can be parameterized in slope/intercept form as: 


exp [Lo=0(Sjo + 7)9i)| 


P(c|6j, dj, 1;) = C1 7 (7) 
Lu=0 EXP [Lo=0(5jo + 7j9i) 
with 6; = [dj0,.. .6¢;-1] essentially representing more than one intercept parameter. 


This divide-by-total model is a constrained version of Bock’s nominal response model 
(e.g., Thissen, Cai, & Bock, 2010), in which there is a single slope parameter for each 
item (or alternatively that the slope parameter for all categories on a single item are 
constrained equal). This model reduces to the 2PL when C; = 2. 

It is easy to see how a monotonic polynomial can be substituted for the linear predic- 
tors of the above item models for a more flexible modeling approach. Liang (2007) was 
the first to describe the logistic function of a monotonic polynomial as simply the 2PL 


model, but with the linear predictor replaced by the monotonic polynomial for item j: 


1 
~~ Dbexp(—1G (6; C0}, T;)) 


P(1|6;,€j, Wj, 0), T;) (8) 

The fact that m;(6;,¢;,j;,«;,T;) is monotonically increasing ensures that the item 
response function for this item is also monotonically increasing. At the lowest-order 
polynomial (i.e., k = 0), this model reduces to the 2PL with the additional constraint 
that the slope parameter is positive. At k > 1, one or more “bends” in the item response 
function may occur, resulting in flexible IRF shapes. 


A generalization of the above model is a GPC model that makes use of the monotonic 
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polynomial in its linear predictor, which we abbreviate GPC-MP, constructed as: 


exp pace + 1; (8, 0), Wj, 7)))| 
Cj-1 


Vueo exp Ea + I (Bj, Cj jy z)))] 


P(c|O;, €j, Wj, 0j, Tj) = (9) 
where m; (Gi, iy ig Ty) 8 DID) ee bongs is the monotonic polynomial without 
the intercept term. The set of C; multinomial logits in Equation (9) share the same 
MP term but have different intercepts. This model reduces to the GPC model when 
k = 0, making m; a linear function of 6. When C; = 2, this model reduces to the 
logistic function of a monotonic polynomial. When both k = 0 and C; = 2, we have a 
2PL model. Furthermore, it can be shown that the conditional response probability of 
category c given the response is in category c or c’ is a logistic function of a monotonic 
polynomial. For the GPC model the intercepts are overparameterized and hence one 
strategy employed by Muraki (1992) for identification is to constrain the first intercept 
Cjo to zero. We take this same approach to identify the GPC-MP model. 

To illustrate the different shapes that the GPC-MP model can take, several example 
IRFs and accompanying parameters appear in Figure 1 and Table 1. These IRFs were 
obtained by fitting models with all items modeled as k = 1,2 or 3 to data later de- 
scribed in Sections 4.1 and 4.3 for dichotomous items (Examples 1-3) and 5-category 
items (Examples 4-6). In general, the order of the polynomial will determine the number 
of possible “bends”’ in the IRFs, though some departures from the 2PL and GPC are 
more noticeable than others. In examining dichotomous items (left column of Figure 1), 
the departures from a 2PL can be mild (Example 1), appear as flat regions in locations 
where we might otherwise expect a lower or upper asymptote before dropping down to 
0 or up to 1 (Example 2), and can be small flat regions or mild kinks within the same 
IRF (Example 3). For polytomous items (right column of Figure 1), wherever there is 
a non-standard bend in the IRF, it will appear in the same region of @ across response 


functions for all categories. Note also that exp(w) is no longer directly interpretable as 
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an overall slope of the item as it is under the 2PL and GPC model, but is interpretable as 
the slope of a tangent line when @ = 0. For instance, Examples 4 and 6 have very small 
values (exp(w) = .17 and .03) where there is a flat region across all IRFs at 9 = 0, but 
Example 5 does not share this property. The remaining a and Tt parameters determine 
curvature changes of the IRF and their interpretation is nontrivial. 
3 Parameter Estimation 

In this section we describe two different estimation approaches for the GPC-MP 
model, beginning with complete data estimation as used in the surrogate-based ap- 
proach used by Liang (2007) and Ramsay (1991). The complete data models are then 
used to derived the EM MML estimation method. Additional details relevant to our 
implementation of the GPC-MP model are also discussed. 
3.1 Complete Data Likelihood and Surrogate-Based Estimation 

Let y;; represent the response from individual 7 to item j. We denote 


Cj-1 
F(yiil@i.nj)) = T] P(cli. 9% (10) 
c=0 


as the conditional density of y;; given the latent trait 0;. All parameters for item j are 
contained in a vector Nis and x-(y;;) is an indicator function that is equal to 1 if and only 
if yij = c, and 0 otherwise. Note that P(c|6;,4;) is a general expression that may be any 
of the item response functions presented in Section 2. Under the local independence 
assumption (Lord & Novick, 1968), the conditional density of person i’s response vector 


y; is then: 
Ff (yil@i, 1) =ILs Yili, 4;)- (11) 


The complete data likelihood is then: 


N 
L(1, #,07|8, Y) = 


If rl) Tle ] (12) 
1 


1= 
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where the bracket on the right-hand side includes the population distribution of 6, which 
in this case is assumed to be normal. For typical single-sample calibration, 1 and 0? are 
fixed for identifiability of the model, yet we retain these parameters in the following 
notation to later facilitate description of multiple group estimation. The log-likelihood 


can be partitioned as 


I(,#,02|8, Y) = log L(]@, Y) + 10g L(y, 07/8) (13) 


where the part pertaining to item parameters corresponds to: 


n N Cj-1 
1(4|0, Y) = log L(q|@,Y) =) | do DY xc(yij) log P(c|6i,4;) | - (14) 
j=l | i=1 c=0 


Note that Equation (14) is a set of n independent multinomial logistic regression log- 
likelihoods. If the latent traits were known we could proceed to estimate item parameters 
by differentiating /(4|@, Y) with respect to each model parameter, and solve the resulting 
system of likelihood equations through iterative methods (e.g., Newton-Raphson). 
Ignoring constants, the part that corresponds to the population distribution mean 


and variance parameters can be written as: 


N N 
I(t, 0°|8) = log L(u,07|8) « — > log(o*) — => Yi(8i— #)”. (15) 


For this complete data model, linear sufficient statistics exist. They are the sum and 
sum-of-squares of the latent trait values: ).., 6; and DN, 6?. 

Prior to estimation of model parameters, Liang (2007) and Ramsay (1991, 2000) com- 
pute an approximation to the latent traits, 86, which are then assumed known in subse- 
quent estimation of model parameters. Although using an altogether different approach 
to estimating IRFs, in Ramsay (1991, 2000) respondents are ranked using weighted 


sum scores, )/j'4 ee. Xc(Yij)Wjc, With weights determined by wjc = logit(PY””)) = 
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logit(P\)), where pe and re represent the proportion of respondents in the up- 
per and lower 25% of the sum score distribution that responded c for item j. Any ties 
among ranks across respondents with different response vectors may be broken ran- 
domly. The ranks are then transformed onto the quantiles of some distribution, usually 
standard normal, and used as provisional estimates of 6. Liang (2007) instead used the 
respondents’ scores on the first principal component (of the covariance matrix of the 
items) as the ranking basis and normalized the resulting scores. Liang (2007) then used 
these scores as approximate @ estimates for use in the optimization of the complete data 
log-likelihood. Following Liang (2007), we refer to this as a surrogate-based (SB) approach. 
3.2 EM MML Estimation 

Our preferred approach is to estimate the GPC-MP model with maximum marginal 
likelihood using the Bock-Aitkin EM algorithm (Bock & Aitkin, 1981). We assume a 
population distribution, such as a normal distribution, $(0;|,07). The marginal (ob- 
served data) distribution of a response vector is obtained by marginalization over the 


unobserved latent traits (Bock & Lieberman, 1970): 


L(y u,7lyi) = flyiln no) = f F(vil6:- 0) (641,07) a0 (16) 


And the marginal log-likelihood is given by I(y, u,07/Y) = a , log L(y, p, o7|y;). 

Integration over 6 must be conducted numerically. For simplicity we employ a Q- 
point rectangular quadrature rule where the quadrature points X1, X2,...Xg are equally 
spaced (e.g., between -6 and 6 on @). Serving as quadrature weights are normalized 
population distribution ordinates, ie., Wy = $(Xq|p,07)/ Eee 19(Xq|u, 07). 

The Bock-Aitkin EM algorithm alternates between two steps. In the E-step, item 
parameters are assumed known and @ is treated as missing data. The conditional ex- 
pectation of the complete data log-likelihood is taken with respect to the posterior 


f(@lyi,4, 4,07). In the M-step, the conditional expected log-likelihood is optimized, 
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resulting in updated parameter estimates. The two steps alternate until only a small 
change in parameter estimates and/or log-likelihoods (e.g., .001) occurs across succes- 
sive iterations. We refer the reader to the Appendices for additional details. 
3.3 A Model with Multiple Groups 

Assuming g = 1,2,...,G groups, the marginal log-likelihood for the combined mul- 


tiple group model is the sum of the within-group marginal log-likelihoods: 
2 2 
L(y, p,0°|Y) = S71 (Hg, Hg, Oe l¥g)- (17) 
gal 


As will be illustrated later in empirical examples, two possible uses of multiple group 
estimation are 1) the ability to fix some (or all) item parameters equal across groups to 
link the scale and estimate whether a studied group’s latent distribution has a different 
mean (and/or variance) versus a reference group with a fixed mean and variance, and 2) 
test for differential item functioning by estimating whether the order of a polynomial for 
a particular item is the same (or different) across groups and whether item parameters 
are equal across groups. Both of these investigations require implementation of equality 
constraints on some item parameters across groups. For this purpose we estimated all 
item parameters simultaneously and used Lagrange multipliers (e.g., Bertsekas, 1996). 

Nested models under this approach can be compared using a likelihood ratio test. 
Similarly, GPC-MP items with higher-order polynomials require the addition of uncon- 
strained « and t parameters and can be compared versus models with lower-order poly- 
nomials using such a test (Liang, 2007). 
3.4 Additional Estimation Details 

In initial tests of the GPC-MP model, some item parameters encountered an ill- 
conditioned likelihood, and thus caused estimation problems and a rank deficient Hes- 
sian. We imposed diffuse priors on a subset of item parameters to aid in estimation 


(e.g., Cai, Yang, & Hansen, 2011; Woods & Thissen, 2006), which technically results in 
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Bayesian estimation of the posterior maximum for model parameters. However, such 
priors may be considered “soft constraints” that serve merely to stabilize estimation as 
is commonly employed with the three-parameter logistic model (Cai et al., 2011); their 
use usually results in trivial changes in the estimated IRF and log-likelihood. Based 
on some preliminary tests of the model, we used diffuse priors of 7(t) ~ N(-—1,500) 
and 7t(«) ~ N(0,500). The choice of a prior mean of 0 for a and a negative value for 
tT roughly correspond to values that would cause the GPC-MP model to reduce to the 
GPC model; although t would need to be —oo for this to occur, small negative values 
appeared to have a similar effect as more extreme negative values. 

All programming and simulations were conducted in R (R Core Team, 2012) by the 
first author. Newton-Raphson with a simple backtracking algorithm was implemented 
for maximization of the complete-data and M-step log-likelihoods: the step size was 
successively halved (until le-6) if the updated parameter estimates did not lead to an 
increase in log-likelihood. Hessians were conditioned by ridging. 

Conceptually, the approach used by Liang (2007) is similar to computing a provisional 
estimate of 6 and then performing a single M-step to estimate the model parameters. It 
should be noted that there were some additional minor differences between our version 
of the surrogate-based approach and Liang’s actual implementation. Specifically, Liang 
did not re-parameterize 6 and A from the monotonic polynomial and instead opted to 
use constrained estimation to enforce all 6B > 0 and A > 0. Furthermore, Liang did 
not utilize diffuse priors on any model parameters and did not ridge the Hessians, but 
instead used the negative information matrix. 

3.5 Model Selection 

Since the GPC-MP model can take a different degree of polynomial, 2k + 1, for each 
item, some selection process may be necessary to avoid blindly picking k. To auto- 
mate this process, we primarily experimented with using Akaike’s Information Crite- 


rion (AIC), which is defined as AIC = —2log L + 2(# of parameters), where log L is the 
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log-likelihood. In the surrogate-based approach, log L is based on the complete-data log- 
likelihood, whereas under EM MML, log L is the marginal log-likelihood. In a subset 
of examples and simulations, we implemented an AIC step-wise approach in which we 
started with a baseline model where all items were k = 0. We looped over all items, 
setting each to k = 1 while keeping all other items the same as in the previous step, cal- 
culated the change in AIC and chose to increase the degree of polynomial for the item 
that improved AIC the most. This process repeated until no additional items modeled as 
k = 1 could improve AIC. The resulting model was used for a starting point in looping 
over all items as k = 2, and so on. For the surrogate-based approach, this step-wise 
model is efficient to compute since 6 estimates are fixed and we would only need to 
loop over all items once; the change in AIC is not dependent upon whether other items 
are modeled with a higher degree polynomial. For EM MML, this process is slower as 
provisional expected @ estimates at each quadrature point may change and modeling 
one item as k = 1, for example, will affect whether modeling a different item to a higher 
degree polynomial will improve AIC. Thus, we also experimented with a more crude 
but efficient approach in which models with all items were estimated as k = 0,1 or 2 in 
simulations (up to k = 3 in empirical examples) and AIC was used to select among these 
models. We refer to these selected models as AIC step-wise and AIC selected, respectively. 
4 Empirical Examples 

The feasibility of the GPC-MP using EM MML is now demonstrated using examples 
from educational and psychological assessments. We show that for empirical datasets 
model fit may improve by including some items as higher-order polynomials. The po- 
tential applications of GPC-MP using EM MML estimation to test mean differences in 
the latent trait across groups and in differential item functioning are also illustrated. 
4.1 Probing for Non-Standard IRFs 

We analyzed responses from 32 items from the Program for International Student As- 


sessment (PISA) Reading Book 8 from the year 2000 using a sample of 3000 participants 
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as described by Cai (2010). Data were in the form of dichotomous scores (incorrect, 
correct) for most items, and partial credit (0,1,2) awarded for three additional items. 
Column 1 of Table 5 shows the number of categories for each item. The purpose of this 
illustration was merely to probe for non-standard IRFs in an educational testing situation 
that included a mix of dichotomous and polytomous items. We estimated five models 
with the first four where all items were modeled as k = 0,1,2,3 and a final model where 
the degree of the polynomial for each item was selected by AIC in a step-wise fashion 
as described in the previous section. 

As shown in Table 2, the best fitting model was the step-wise model, followed by the 
k = 1 model. The GPC-MP model resulted in modest improvement over standard GPC 
and 2PL (i.e., k = 0) items, and over half of the AIC step-wise items were modeled with 
higher order polynomials: 8 items modeled as 5th-order polynomials and 13 items as 
3rd-order polynomials (see Table 5). Estimated item response functions under the AIC 
step-wise model appear in Figure 2. Some k = 1 items had merely a small kink in the 
IRF (e.g., Items 2, 7, 26, 30, and 32) whereas k = 2 modeled items had greater non-trivial 
departures from a logistic IRF (e.g., Items 1, 4, 6, 9, 15, 20, 23, and 24). 
4.2 Estimation of Group Mean Differences 

The calibration data from the Patient Reported Outcomes Measurement Information 
System (PROMIS®) smoking module (Hansen et al., in press; Shadel, Edelen, & Tucker, 
2011) included responses from 4,201 daily smokers (28-30 days smoked out of the past 
30 days) and 1,183 non-daily smokers (<28 days smoked out of the past 30 days) to items 
divided among 6 different constructs related to smoking behavior: nicotine dependence, 
hedonic benefits, coping benefits, social benefits, psychosocial risks, and health risks. 
Data collection followed a randomized block design using 26 overlapping forms to re- 
duce respondent burden. Items had five possible ordered response options. 

For illustration purposes, we analyzed responses to only the hedonic benefits item 


bank, which included 16 items for daily smokers and 17 items for non-daily smokers, 
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15 of which were common to both groups. Thus, the goals of this illustration were to 
1) Model the 15 common items using the GPC-MP model and AIC to select the order 
of the polynomial for each item; and 2) Estimate the mean difference in daily /non-daily 
smokers on the underlying hedonic benefits latent trait. 

To accomplish these goals, the AIC step-wise approach described in Section 3 up to 
k = 3 was used while the daily smoker mean and variance of hedonic benefits were 
fixed to a standard normal distribution and the non-daily mean and variance was freely 
estimated. Since the original calibration analyses suggested that the 15 common items 
differed negligibly in their item parameters across groups, parameters for these items 
were constrained equal across groups. The items not shared by both groups were fitted 
with the GPC model. 

As shown in Table 5 and Figure 3, AIC selection of k for each item suggested several 
departures from the standard GPC model, including seven items with a 7th-order poly- 
nomial (k = 3). The final model also indicated that non-daily smokers had a lower mean 
on hedonic benefits (fi = —.42) than the daily smokers (fi = 0), x2(1) = 76.72, p < .001 
(see also Table 3). 

4.3 Differential Item Functioning 

The GPC-MP model using EM MML estimation also may have utility in situations 
where the order of polynomial for a particular item may differ across groups, and may 
additionally be used to test whether higher order polynomials have the same item pa- 
rameters across groups. For this illustration, we used the social benefits domain from 
the PROMIS® smoking module, which included 12 items for daily and non-daily smok- 
ers, 9 of which were common to both group’s item bank. Although initial calibration 
results indicated few items with substantial differential item functioning as indicated by 
comparing the expected score of the IRFs across groups (Hansen et al., in press), for 
illustration purposes we selected only two anchors (items 5 and 7) modeled as standard 


GPC items. The remaining 7 items were studied one-by-one using the GPC-MP model 
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while freely estimating the mean and variance of the reference group’s latent trait. 

One possible procedure is to first determine the order of polynomial for each group 
by fitting a series of models whereby the studied item is k = 0,1,2,3 for a single group, 
and use likelihood ratio tests to determine whether any of the higher order polynomials 
presented significant improvement over lower-order polynomials.* Although this does 
not explore all possible combinations of k for the item across groups, as with use of AIC 
this process should provide a reasonable guess as the to order of polynomial for each 
group as supported by the data. 

The above procedure resulted in only item 1 being modeled as the same order poly- 
nomial (k = 1) across groups (see also Table 5). Constraining all 7 parameters for this 
item equal across groups resulted in significant worsening of model fit, x7(7) = 30.30, 
p < .001, suggesting differential item functioning for this item. Constraining only the 
slope and bend parameters, w, «, and Tt equal across groups did not result in worse fit, 
x7(3) = 2.13, p = .55, analogous to failing to provide evidence for non-uniform differ- 
ential item functioning. Additionally constraining intercept parameters, ¢, equal across 
groups resulted in significant worsening of model fit, x7(4) = 28.17, p < .001, analogous 
to providing evidence for uniform differential item functioning. 

5 Simulation Studies 

Two simulation studies were performed. The goal of the first was to test item re- 
sponse function recovery of the GPC-MP model using EM MML under a variety of 
conditions and against the standard 2PL and GPC model. The second simulation was 
intended as a more focused comparison of GPC-MP using EM MML versus the alterna- 
tive IRF estimation procedures introduced by Liang (2007) and Ramsay (1991). 

The main outcome of these simulations was IRF recovery as measured by a general- 


ization of Root Integrated Mean Square Error (RIMSE; e.g., Liang, 2007; Ramsay, 1991) 


We implemented this approach using a conventional p < .05 threshold for hypothesis testing though 
note that Benjamini-Hochberg adjustment is sometimes used in practice to control the false discovery rate 
in differential item functioning situations (e.g., see Thissen, Steinberg, & Kuang, 2002). 
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for polytomous item j: 


Fo (£5)(8q) — ES;(0,))(8q) \ 


Se (9) 


RIMSE; = x 100 (18) 


where ES;(6,) = Ey c- P(c|6q,4;)/(Cj — 1) is the expected score for item j at a given 
level of @ (divided by C; — 1 to place dichotomous and polytomous items on the same 
scale), the sum is across quadrature points, and ¢ is the density function for a standard 
normal distribution. The secondary outcome was latent trait recovery as measured by a 


similar index: 


N 1/2 
RMSE¢ = ( y (6; — 0?) x 100 (19) 
i=1 


where 6; is the actual latent trait for subject i and 6; is the estimated trait.? Unless noted 
otherwise, trait estimates were calculated using the expected a posteriori (EAP) scoring 
method (Bock & Mislevy, 1982). Conceptually, RIMSE; measures the agreement of the 
estimated IRF with that of the true IRF with more weight given to estimates near the 
center of the latent trait distribution. To provide a single index of IRF (and trait) recovery 
for each cell in the simulation design, RIMSE; (and RMSEg¢) was further aggregated 
across items (or respondents) for each replication and again across all replications in a 
single cell. 
5.1 Study 1: IRF Recovery Versus Standard Approaches 
5.1.1 Data Generation 

In assessing the GPC-MP model’s IRF recovery, we generated unidimensional data 
crossing the following data generation conditions: 1) Number of categories per item (2 
and 5), 2) True polynomial order (k = 1 and k = 2), 3) Number of test items (10 and 
20) and 4) Sample size (500 and 3000). One-hundred data sets per cell were generated. 


In all cases, the latent traits were drawn from a standard normal distribution. For half 


3We thank an anonymous reviewer for suggesting this index for trait recovery. Results for latent trait 
recovery using RIMSE, as in Liang (2007) are available from the authors upon request. 
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of the test items, the true item parameters were generated randomly across replications. 
For dichotomous items, € ~ unif(—1,1), exp(w) ~ unif(.3,2.5), « ~ unif(—1,1), and 
exp(T) ~ unif(0,1), which are similar to values used in simulations by Liang (2007). 
For polytomous (5 category) items, > = 0, ¢, ~ unif(1,2.5), ¢2 ~ unif(—1,1), é3 ~ 
unif(—1,1), 4 ~ unif(—2.5,—1) and exp(w) ~ unif(.3,1.5). 

The true parameters for the remaining items were the same across replications. Di- 
chotomous item parameters for k = 1 and k = 2 were drawn from initial calibration 
results for five different items estimated using a subset of the PISA dataset. Polytomous 
item parameters were drawn from five items of the PROMIS® social benefits domain as 
calibrated on daily smokers only. For the 20 item cell, parameters from five items were 
repeated to obtain parameters for 10 items. Selection of these item parameters was based 
on which items had k = 1 or k = 2 models that improved AIC and a visual inspection to 
include extreme-looking items and a variety of IRF shapes. 

5.1.2 Fitted Models 

In each cell of the design, EM MML was used to fit a series of GPC-MP models. This 
included three models that modeled all items using the same degree polynomial from 
k = 0,1,2, and a fourth AIC selected model that merely picked among the all k = 0,1 
or 2 models. In a subset of conditions (10 items), we included use of the AIC step-wise 
model starting at k = 0 and working up to k = 2. All models were estimated using 
49 equally spaced quadrature points between -5 and 5, and estimation terminated upon 
changes in item parameters less than .001. 

5.1.3 Results 

RIMSE results for items appear in Figure 4; lower RIMSE values are better. To facil- 
itate comparisons among fitted models, for any given data generation combination, the 
fitted models were ranked (1 = best, 5 = worst) with ties allowed and shading assigned 
to these ranks (white = better, gray = worse). As our key interest is whether higher- 


order polynomials have better IRF recovery than lower-order polynomials and whether 
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AIC selection performed well, this shading strategy emphasizes differences across fitted 
models, but not across data generation conditions. 

Even when the true model was a higher-order polynomial (e.g., k = 2), IRF recovery 
was not always best by fitting the true model to every item, but this trend depended 
on the number of categories per item. For dichotomous items, IRF recovery was often 
better with a lower-order polynomial (e.g., k = 0 or k = 1). This was especially the case 
with few respondents and items. For instance, IRF recovery at N = 500 and 10 items 
was best with k = 0 when the true k = 1 (RIMSE = 5.14) and k = 1 when the true k = 2 
(RIMSE = 5.43). At N = 3000 and 20 items, however, IRF recovery was best when the 
true model was fit to the data (RIMSE = 1.75 at k = 1 and RIMSE = 2.59 at k = 2). 
With polytomous items, IRF recovery appeared to be at or near best when fitting the 
true model, regardless of number of items and sample size. For example, IRF recovery 
was always first or second best among fitted models when all k = 1 items were modeled 
as k = 1, and the same for when k = 2 items were modeled as such. 

With respect to the utility of AIC in selecting the degree of each items’ polynomial, 
the step-wise model tended to be third best (6/8 cells) or better in the conditions in 
which it appeared. Although this performance appears modest at best, note also that the 
AIC step-wise model never had RIMSE values more than .3 behind the top performing 
model, whereas the worst-performing models could have RIMSE values that were 3 or 
more behind. With the exception of only two cells in the design in which it was less 
than .1 behind the second place model (both at N = 500 and 10 items), the cruder AIC 
selected model was always either first or second in IRF recovery. The AIC selected model 
performed exceptionally well in conditions with more items and respondents, in which 
it was always first or tied for first with RIMSEs ranging from .96 to 2.59. 

Recovery of latent traits as measured by RMSEg nearly mirrored recovery of IRFs (see 
Figure 5) in that latent trait recovery tended to be best in cases when IRF recovery was 


best. For example, under 10 items and N = 3000, the ranking of the methods was iden- 
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tical to IRF recovery for all cells but three. Although the ranking of methods departed 
slightly more from IRF results under other conditions, it was still true that latent trait 
recovery tended to be better for higher-order polynomial models under conditions of 
more information (polytomous items, more items, more subjects). Results for the AIC 
selected model were similar in that it was ranked first or second in all cases but one at 
N = 500 and 10 items. The AIC step-wise model also performed decently in terms of 
latent trait recovery, being ranked third in 5/8 cases, and first under two other condi- 
tions. Overall these results suggest that use of the GPC-MP model may be best with a 
high degree of information, but also that the use of AIC can aid in choosing the degree 
of the polynomial for each item when the true data generating model is not known. 

5.2 Study 2: Comparison Versus Alternative Estimation Approaches 

5.2.1 Data Generation 

In order to compare the GPC-MP model to alternatives, we crossed the following data 
generation conditions: 1) Number of categories per item (2 and 3), and 2) Sample size 
(500 and 3000). In each cell of the data generation design, one-hundred datasets with 
twenty items each were generated. Data generation differed from that of the previous 
study in that the true IRFs were not from GPC-MP items, but were instead generated 
from a mixture of the cumulative distribution function (CDF) of normal variables so that 
our approach is not implicitly favored. 

Half of the IRFs under dichotomous items were generated following Liang (2007) 
who used two normal CDFs, piN (11,07) + poN(p2, Ga); with pi ~ unif(.3,.7), po = 
1— py, py ~ N(-1.5,.17), 2 ~ N(1,.17), 01 ~ N(1,.12), 02 ~ N(.4,.17). In our subjective 
opinion, the resulting IRFs under such data generating values were often very similar 
and not very extreme. Thus, the remaining dichotomous items were a mixture of three 
CDFs with p; and pz ~ unif(.1,.4), p3 =1— pi — pr, fi ~ N(—1.5, 47), fo ~ N(1.5,.4), 
w3 ~ N(0,.47), and 01,02, and 03 independently drawn from N(.4,.17). Polytomous 


(3 category) items were generated using a mix of 2 or 3 (half of the items for each) 
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normal CDFs pieced together in an analogous fashion to a graded response model, i.e., 
Ply 2c) = piN(pat te, oz) + poN (2 + te, 05), with P(yj = 0) =1and 4, > t¢_1 since 
ty ~ unif(—1.5,—.5) and i) ~ unif(.5,1.5). All # and o values were drawn using the 
same distributions as dichotomous items. 

5.2.2 Fitted Models 

The main goal was to compare the following approaches: 1) GPC-MP using EM 
MML,; 2) GPC-MP using the surrogate-based approach; and 3) The kernel smoothing 
approach of Ramsay (1991). To this end, GPC-MP models were estimated in the same 
manner as in the previous simulation study (all k = 0,1,2 models, the AIC step-wise 
model, and AIC selection of k = 0,1,2 models). In addition, these models were estimated 
using EM MML and the surrogate-based approach. Finally, the ksIRT function from 
the KernSmoothIRT (Mazza et al., 2013) package in R was used to estimate IRFs under 
the kernel smoothing approach using 51 equally spaced points between -3 and 3 on @. 
Note also that ksIRT generates only maximum likelihood-based latent trait estimates; we 
report the ML estimates but also added a normal prior for computation of EAP estimates 
for this method. 

Theoretical expectations about the expected pattern of results are difficult to form. 
Note that our implementation of the EM and SB approaches are the same except that the 
latent traits, 0, are treated as missing data under EM but approximate estimates are used 
in the complete data likelihood under SB. To our knowledge, the performance of such 
initial latent trait estimates under SB has not been studied versus the EM MML solution. 
Had further iterations of the SB approach been employed alternating between trait and 
item parameter estimates, it may resemble the joint maximum likelihood approach. The 
KS approach also uses an altogether different approach to initial 6 estimates and a non- 
parametric approach to estimating IRFs. Finally, the data generating IRFs are not derived 
from any of the estimation approaches we employed, and the use of a graded-type 


(i.e., differencing) model to piece together IRFs for polytomous items may present an 


October 30, 2014 23 


additional challenge for the GPC-MP items (which are divide-by-total models). 
5.2.3 Results 

IRF recovery for the various modeling approaches appears in Figure 6 using the 
same shading strategy as in Study 1. At N = 500, the KS approach emerged as the 
top performing method by .08 and .23 over the second best method (under dichotomous 
and polytomous items, respectively), with the remaining results somewhat mixed. In 
examining all k = 0,1,2 for both SB and EM MML, the all k = 1 models performed well 
under dichotomous items at that sample size, suggesting that higher-order polynomials 
were appropriate for the IRFs, but the k = 0 model performed best under polytomous 
items. In comparing the utility of AIC, for dichotomous items, EM MML (k = 1 and AIC 
selected) and SB with the AIC selected model were tied for second with dichotomous 
items, followed by the AIC step-wise model for EM MML. For polytomous items, AIC 
step-wise and selected for EM MML performed second and fourth best, respectively, 
suggesting that EM MML combined with some form of AIC selection is viable. At 
N = 3000, the advantage of higher-order polynomials and EM MML over SB emerged 
more clearly. With all k = 2, EM MML occupied the first and second place (tied) spots 
for dichotomous and polytomous items, respectively, with the EM MML AIC step-wise 
and selected models occupying places between second and fourth. The best competing 
method was the KS approach, which was .72 worse than the EM step-wise model under 
dichotomous items, but .11 better under polytomous items. 

Recovery of latent traits presented a similar pattern to that of IRF recovery with a 
few exceptions (see Figure 7). At N = 500, KS with EAP scoring performed better 
by .02 under dichotomous and by 2.49 under polytomous items over the second best 
method. Under dichotomous items at N = 500, SB and EM MML with higher order 
polynomials (anything but k = 0) or AIC step-wise/selection performed well, but with 
SB having a slight advantage. With polytomous items at N = 500, k = 0 models per- 


formed better than those with higher-order polynomials and AIC selection/step-wise 
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models performed better for EM MML than for SB. As with IRF recovery, at N = 3000 
the advantage of the GPC-MP model with EM MML began to emerge along with a 
clearer pattern of results. The k = 2, AIC selected, and AIC step-wise models for EM 
MML occupied the top two spots (with ties) along with SB with AIC selection. Under 
polytomous items, these same models occupied spots behind the KS method with EAP 
and ML scoring. ML estimates using the KS approach, the default estimates provided 
by KernSmoothIRT (Mazza et al., 2013), tended to perform worse than other methods 
and was always ranked last with the exception of its second place performance under 
N = 3000 with polytomous items. 

6 Discussion 

We have presented a monotonic polynomial generalized partial credit model that 
subsumes the 2PL and GPC models and is estimated using EM MML. This model has the 
potential to probe for and model non-standard item response functions and has potential 
applications in multi-group analyses. In addition, our simulation studies demonstrate 
that the use of higher order polynomials can result in better IRF and latent trait recovery 
when the true IRF does not follow a 2PL or GPC shape. Furthermore, the GPC-MP item 
model using EM MML estimation tended to perform better than the SB approach of 
Liang (2007) for IRF recovery and was nearly on par with the IRF recovery of the Kernel 
smoothing approach of Ramsay (1991). 

That said, we note some important limitations and potential directions for future 
research. First, as shown in Study 1, use of a lower-order polynomial can sometimes 
lead to better IRF recovery than use of a higher-order polynomial even when the true 
IRF follows the higher-order polynomial. We speculate that when insufficient informa- 
tion exists (e.g., few items, subjects), that the use of higher-order polynomials is prone 
to fitting noise in a similar way to how the kernel smoothing approach may result in 
non-monotonic IRFs (especially in the tails) even when the true IRF is monotonically 


increasing. Although not explicitly manipulated in our studies, we would also suspect 
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that more discriminating items may aid in providing more information necessary for 
fitting the GPC-MP. 

To enhance the usability of the GPC-MP, a mechanism may be required for either 
probing for higher-order polynomials and/or testing whether a particular item is better 
modeled as a higher-order polynomial. In our simulations, AIC demonstrated promise 
in selecting the order of polynomial. AIC selected and step-wise models were not al- 
ways the best out of fitted models, but were never far behind and are clearly preferable 
to blindly picking the order of the polynomial. Although our choice of AIC is somewhat 
arbitrary, Liang (2007) experimented with use of AIC, BIC and likelihood ratio tests, find- 
ing that these approaches performed similarly. We suggest that any of these approaches 
remain as a potential direction for additional research, and add that examination of item 
fit statistics (Orlando & Thissen, 2000) also remains a possibility. We speculate that BIC 
would result in fewer items being modeled as higher-order polynomials due to its bias 
towards parsimony, and would thus allow researchers to model more items as 2PL or 
GPC items and reduce computational time if a BIC step-wise model were implemented. 

Under our simulations, better IRF recovery for the GPC-MP model typically trans- 
lated into better latent trait recovery. One exception to this pattern was at N = 500 with 
dichotomous items for Study 2, where better IRF recovery using EM MML did not nec- 
essarily lead to better trait recovery than the SB approach. Since these two approaches 
both used the GPC-MP model, without additional research it is difficult to determine 
whether there is a bias/efficiency trade-off to latent trait recovery under these two esti- 
mation approaches or if another explanation exists. To our knowledge, the SB approach 
used by Liang (2007) has not been extensively studied in previous research in compari- 
son to EM MML, nor has the use of provisional @ estimates by Ramsay (1991). We note 
that the current studies were primarily geared towards a comparison of IRF recovery 
and not specifically designed with an investigation of trait recovery in mind. Thus, as 


the KS approach tended to perform near the top in terms of IRF recovery under Study 
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2, we also do not know if this is due to kernel smoothing of IRFs, surrogate @ estimates 
that are better than Liang’s use of the subjects’ scores on the first principal component, 
or whether piecing together CDFs as in the graded response model instead of a divide- 
by-total manner put the GPC-MP at a disadvantage for polytomous items. Nonetheless, 
even though true items were not generated by the GPC-MP model, the GPC-MP model 
using AIC step-wise or selection was never far behind the KS approach, and sometimes 
outperformed it. 

Several possible variations of the GPC-MP model and comparisons with other ap- 
proaches also remain for future research. For instance, should researchers wish to es- 
timate a GPC function of a polynomial that has a negative relationship with the latent 
trait or without the monotonicity requirement, constraints on A or both A and 6 could 
be released, respectively, by not reparameterizing. While we do not rule out the possi- 
bility of constructing GPC functions of other types of polynomials (e.g., splines), such 
an approach will likely not offer the same level of parsimony as GPC-MP, and it is un- 
clear whether this would have any advantage over existing approaches (e.g., monotone 
splines; Ramsay & Winsberg, 1991). Finally, our approach has not yet been compared 
with other promising, though computationally intensive Bayesian non-parametric ap- 
proaches (e.g., Duncan & MacEachern, 2008; Miyazaki & Hoshino, 2009).4 

In conclusion, use of the GPC-MP using EM MML has potential advantages over the 
KS and SB approaches. In the realm of multiple-group analysis, EM MML provides a 
more natural way of estimating group differences and testing of differential item func- 
tioning across groups. Whereas the KS approach models all items non-parametrically, 
the GPC-MP model has the potential to allow researchers to model some or most items 
on a test as standard 2PL and GPC items. Thus, use of the GPC-MP allows greater par- 
simony, potentially easier interpretation of item parameters in a test, and formal testing 


of whether a higher-order polynomial fits better than a lower-order polynomial. 


4See Liang (2007) for a comparison of a logistic function of the monotonic polynomial using surrogate- 
based estimation with the methods of Duncan and MacEachern (2008) and Qin (1998). 
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Appendix A. Further Derivations of EM MML Estimation 
Recall that the complete data likelihood is L(y, p,07|y;,0;) = f(yi{@:,7)0(0;|1,07) 
for individual i. Thus, an individual’s contribution to the marginal likelihood can be 


approximated to arbitrary precision as 
: »n € 
Li(q,u,o") = Df (yilXq. 9) Wo. 
qa 


Using quadrature, the height of the complete data likelihood at quadrature point Xj can 
be represented as 


n G1 
j=1 c=0 
This suggests that the ordinate of the posterior f(@lyi,, 1,07) at quadrature point X, 


can be approximated by 
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For the item parameter part, the conditional expectation of log L(7|6;, y;) is 
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where 1/,, j/x, and 72 are the current/provisional parameter estimates. Using Equation 


(20), this conditional expectation may be approximated by quadrature as 
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Summing over all N individuals and rearranging terms in the summation, the condi- 
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tional expectation of log L(y|6, Y) from Equation (14) is 
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where Figo = eae Xc(yij) Pi(Xq) is the conditional expected frequencies for item j, cate- 
gory c, at quadrature point q. Item parameter estimates are updated in the M-step by 
treating the expected frequencies as weights and maximizing E(y|1,, x, 02). 

For distributional parameters for the latent traits, pp and 0%, if not fixed in a multiple 
group analysis, can be also be estimated by calculating the mean and variance of the 
expected counts (e.g., see Baker & Kim, 2004). This is possible because the conditional 
expectations of the linear sufficient statistics may also be approximated via quadrature 


and the M-step is closed-form. For example, the conditional expectation of Y™, 6; is 
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where fy = ae P;(X,) is the conditional expected frequencies at quadrature point q. 
An updated latent variable mean estimate is therefore N~! eg 7,Xq. For the variance 


parameter, the conditional expectation of DN , 6? is 
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An updated variance estimate is found by N~! Bele 7X; 3 (N = Eee 7X) . 
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Appendix B. Complete Data Derivatives for the GPC-MP Model 
The complete data log-likelihood for a single GPC-MP item ] is: 


N 
= y ye Xe(yis) log P(c|6j, ¢j, Wj, &), Tj): 


where P(c|6j,¢;,;,«;,T;) is the item response function for item j under the GPC-MP 


model. Differentiating with respect to a typical parameter, 77;, leads to: 
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where mj; = = mj (6, W;,@;,T;) is short-hand for the monotonic polynomial without the 


J 

intercept parameter, and Pj, = P(ul6;, Gir Wi, Oj, T;) is short-hand for person i’s probabil- 
am, 

ity of ans to category u on item j under the GPC-MP model. Of course, a =0, 

and Se is 1 when differentiating with respect to ¢j, and 0 otherwise. The derivatives 


of mi, i the parameters (vj, aj, and Tj, are simply the following and can be substituted 


into the above equation: 
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in which s = 1,2,...,k, the T matrices are specific to item j, and -5" is the vector: 
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The derivatives of the matrices T with respect to aj; and Tj; have the form: 


0 0 0 0 0 0 
—2 0 0 0 0 0 
2aj, —2 0 6° or 66 
9 0 2a, —2 O° oO O 
djs s _ : , 
0 90 © ass Oa, =2 0 
0 0 0 0 2w, —2 
0 0 0 0 0 2a, 
0 0 0 0 0 0 
0 0 0 0 0 0 
exp (Ts) 0 0 0 0 0 
ca 0 exp(Tjs) 0 0 0 0 
s= 
OT;s . 
0 0 0 exp(Tjs) 0 0 
0 0 0 -:- 0 exp(Tjs) 0 
0 0 0 ::- 0 0 exp(Ts) 


For computing the Hessian, we used the following cross-product of gradients approx- 


imation, in which derivative vectors are computed for each individual and their outer- 


oe oi; 
The above derivatives can be easily adapted for the M-step computations in EM MML 


/ 
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products are accumulated (e.g., Bock & Lieberman, 1970): DN, ( AU ‘) ( ds ") 


estimation by summing over the Q quadrature points instead of over the N individuals, 


and by treating the expected frequencies 7j,- at each quadrature point as weights. 
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Table 1: Item Parameters from Example IRFs 


Example 1 Example 2 Example3 Example4 Example5 Example 6 


C1 
Go 
C3 
C4 
Ww 
ay 
x2 
a3 
al 
TM 
T 


0.37 0.71 -0.90 0.14 3.48 -0.21 
0.18 2.50 -0.24 

-1.06 -0.25 -1.20 

-1.52 -1.64 -0.93 

-0.11 0.69 0.54 -1.76 1.02 -3.41 
0.24 -0.50 -0.73 3.80 0.89 -0.53 
0.52 0.81 -0.74 0.52 

0.36 8.20 

-0.21 -8.48 -6.65 -1.79 -8.70 -9.92 
-3.32 -1.96 -8.99 -5.91 

-8.26 -1.50 


Note. Examples 1 to 3 are dichotomous items with k = 1,2 and 3, respectively. Examples 
4 to 6 are 5-category items with k = 1,2 and 3, respectively. 
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Table 2: Model Overview for PISA data 


Allk=0(GPC) AllLk=1 Allk=2 Allk=3 AIC step-wise 


36 


# Parameters 67 131 195 259 125 
—2logL 106007.9 105761.1 105644.8 105597.6 105704.6 


AIC 106141.9 106023.1 106034.8 106115.6 105954.6 


Note. Best two AIC values are underlined. 
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Table 3: Model Overview for PROMIS® Hedonic Benefits data 


Allk =0 (GPC) AIC step-wise Constrained Means 


# Parameters 92 158 157 
—2logL 124715.3 124296.2 124373.0 


AIC 124899.3 124612.2 124687.0 


Note. In all models common item parameters were constrained equal across daily and 
non-daily smokers. Group means were constrained equal in the “constrained means” 
model, but the non-daily mean was free in the other models. 
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Table 4: DIF models studying item 1 for PROMIS® Social Benefits data 


Unconstrained Equalw,a,andt Equalw,a«a,t, and ¢ 


# Parameters 116 113 109 


—2logL 88640.2 88642.4 88670.5 


AIC 88872.2 88868.4 88888.5 


38 


October 30, 2014 39 


Table 5: Item Overview for PISA and PROMIS® data 


PISA Hedonic Benefits Social Benefits 

Item Cj; Selected k Cj; Selectedk C; Selected k; Selected ky 
1 2 5 2 5 1 1 
2 5 2 5 1 0 
3 5 1 5 1 0 
4 5 1 5 0 1 
5 5 1 5 0 0 
6 5 3 5 1 2 
7 5 1 5 0 0 
8 5 3 5 1 0 
9 5 3 5 0 1 

5 3 

5 3 

5 3 

5 3 

5 2 

5 2 


NOrPP RPP PRP PPE 
SOON DVT KFPWNF O 


NNNNNNN N 
CON DU FPWN FH 


WN 
© oO 
NNNNNNNNNNNNNNNNNNNWWNNNNWNNNNNN 


Soe eS RK Oe ON NR Re io HK Oo ON OOF oO oN Oe NRE Ne 


Wo 
Nem 


Note. For PROMIS® data (hedonic benefits and social benefits), only items common to 
both daily and non-daily smokers are presented. k; refers to daily smokers and kp refers 
to non-daily smokers. Items 5 and 7 for social benefits served as anchors. C; is number 
of categories per item. k was selected using AIC for PISA and hedonic benefits, and 
likelihood ratio tests for social benefits. 
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Figure 1: Example IRFs 
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Note. Response curves for higher categories are darker. 
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Figure 2: IRFs for PISA AIC Stepwise Model 
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Figure 3: IRFs for PROMIS® Hedonic Benefits AIC Stepwise Model 
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Note. Response curves for higher categories are darker. 
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Figure 4: RIMSE for Study 1 item response functions 


10 Items 
, 5.0 
All k=0 45 
All k=1 5.43 40 
All k=2 
@ AIC Selected 5.15 5.71 3.5 
2 AIC Step-wise SAlié 5.62 3.0 
3 All k=0 oh 
in All k=1 
All k=2 2.0 
AIC Selected 1.5 
AIC Step—wise 
1.0 
k=1 k=2 k=1 k=2 
Data Generating Model 
20 Items 
4.0 
All k=0 3.5 
All k=1 
a All k=2 3.0 
oO 
8 AIC Selected oe 
3 All k=0 
F= 2.0 
LL All k=1 
All k=2 1.5 
AIC Selected 
1.0 


k=1 k=2 k=1 k=2 
Data Generating Model 


Note. Lower RIMSE values are better. Shading (white = better) indicates rank (1 = best) 
of each fitted model’s performance for a particular data generation combination versus 
other fitted models. 
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Figure 5: RMSEg for Study 1 latent trait estimates 


10 Items 
; 5.0 
ako | | 8766 48 
All k=1 48.65 47.28 a 
All k=2 
@ AIC Selected 48.67 47.44 3.5 
2 AIC Step—wise 48.68 47.49 3.0 
. All k=0 ae 
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All k=2 4] 4732) 45.73 a0 
AIC Selected 46.95 45.69 15 
AIC Step-wise 47.10 ne 
k=1 k=2 k=1 k=2 
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All k=1 
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® All k=0 
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k=1 k=2 k=1 k=2 
Data Generating Model 


Note. Lower RMSE values are better. Shading (white = better) indicates rank (1 = best) 
of each fitted model’s performance for a particular data generation combination versus 
other fitted models. 
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Figure 6: RIMSE for Study 2 item response functions 
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Note. Lower RIMSE values are better. Shading (white = better) indicates rank (1 = best) 
of each fitted model’s performance for a particular data generation combination versus 
other fitted models. KS = Kernel smoothing; SB = Surrogate-based estimation; EM = EM 
MML (Maximum marginal likelihood using the EM algorithm). 
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Figure 7: RMSEg for Study 2 latent trait estimates 
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Note. Lower RMSE values are better. Shading (white = better) indicates rank (1 = best) 
of each fitted model’s performance for a particular data generation combination versus 
other fitted models. KS = Kernel smoothing; SB = Surrogate-based estimation; EM = EM 
MML (Maximum marginal likelihood using the EM algorithm). 


